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Bubble  or  drop  distortion  in  a  straining  flow  in  two 
dimensions 
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Joseph  B.  Keller 

Departments  of  Mathematics  and  Mechanical  Engineering,  Stanford  University,  Stanford, 
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(Received  21  January  1980;  accepted  16  May  1980) 

The  distortion  of  a  two-dimensional  bubble  (or  drop)  in  a  straining  flow  of  an  inviscid  imcompressible 
fluid  is  examined  theoretically.  Far  from  the  bubble  the  stream  function  of  the  flow  is  assumed  to  be  axy , 
where  a  is  a  constant.  Within  the  bubble  the  pressure  is  assumed  to  be  a  constant  pb,  and  the  bubble 
surface  is  assumed  to  have  a  surface  tension  a.  Then,  the  shape  of  the  bubble  depends  upon  the  single 
dimensionless  constant  y  =  2(pb  -  ps)/(laa)2lipu\  where  p  is  the  fluid  density,  pt  is  the  stagnation 
pressure  of  the  flow,  and  the  size  of  the  bubble  is  proportional  to  (2 a/pa2)x,i.  For  y  large,  it  is  found 
that  the  bubble  tends  to  a  circle  of  radius  (2<7//kz2)W3  y“‘.  As  y  decreases,  numerical  solutions  show  that 
the  bubble  at  first  becomes  a  square  with  rounded  comers.  Then,  it  develops  four  horns  or  spikes  with 

large  curvature  near  their  ends.  Finally  at  y - 1.8,  the  two  sides  of  each  spike  touch  each  other  near 

the  tip  and  enclose  a  small  bubble  there.  It  is  also  found  that  there  is  a  maximum  value  of  the  Weber 
number  above  which  there  is  no  steady  solution. 


I.  INTRODUCTION  AND  FORMULATION 

In  the  mixing  of  two  fluids,  a  drop  or  bubble  of  one 
fluid  will  be  distorted  and  possibly  split  into  smaller 
parts  because  of  the  flow  of  the  other  fluid  around  it. 

In  order  to  study  this  phenomenon,  we  consider  the 
two-dimensional  case  of  a  drop  or  bubble  of  one  fluid 
in  a  steady  flow  of  another  fluid,  assumed  to  be  inviscid 
and  incompressible.  We  take  into  account  the  surface 
tension  <7  at  the  interface,  but  we  ignore  the  flow  inside 
the  drop  or  bubble,  assuming  that  the  pressure  is  a 
constant  f)b  throughout  it.  From  now  on  we  shall  write 
“bubble"  to  mean  either  bubble  or  drop. 

In  order  to  formulate  this  problem  we  assume  that  the 
stream  function  of  the  flow  without  the  bubble  is  ary, 


finding  the  flow  consists  of  determining  x  +  iy  as  an 
analytic  function  of  (p  +  tip  in  the  half-plane  ipz  0  satis¬ 
fying  Eq.  (1)  at  infinity.  Then,  the  bubble  surface  is 
given  by  setting  ip  =  0  in  x(<p  +  i0)  and  y((p  +  iip),  and  letting 
<p  range  from  -  £  to  The  symmetry  conditions  re¬ 
quire  that  the  bubble  surface  be  normal  to  the  x  and  y 
axes  where  it  intersects  them,  which  yields 

n(~£>0)=r„(i0)=0  .  (2) 

On  the  bubble  surface  the  pressure  in  the  fluid,  which 
is  given  by  the  Bernoulli  equation,  must  differ  from  pb 
by  uk ,  where  k  is  the  curvature  of  the  interface.  This 
leads  to  the  boundary  condition 


where  a  is  a  constant  and  x  and  yare  Cartesian  coor¬ 
dinates.  This  flow  is  symmetric  about  both  the  x  and  y 
axes,  so  we  assume  that  the  bubble  and  the  flow  around 
it  will  have  the  same  symmetry.  Then,  it  suffices  to 
determine  the  flow  only  in  the  first  quadrant  x  s*  0, 
y>  0  (see  Fig.  1). 

We  introduce  dimensionless  variables  by  choosing 
(20/pa2)1 /3  as  the  unit  length  and  (2cra/p)l/3  as  the  unit 
velocity.  We  also  introduce  the  dimensionless  potential 
function  b<p  and  stream  function  6^.  Here,  6>0  is  a 
dimensionless  constant  to  be  chosen  so  that  <p  =  5  and 
<p  =-  \  at  the  stagnation  points  on  the  x  and  y  axes,  re¬ 
spectively.  We  denote  the  streamline  along  the  two 
axes  and  along  the  arc  of  the  bubble  boundary  in  the 
first  quadrant  by  0  =  0.  In  these  variables  bip~xy  and 
b<p  ~(x?  -y2)/ 2  at  infinity,  so  that  +  *y)a/2 

or,  equivalently, 

x  +  iy  ~(2 b)1/2(tp  +  (1) 


at  infinity.  The  flow  occupies  the  region  0  2  0  of  the 
<P,  0  plane,  and  the  bubble  boundary  corresponds  to  the 
segment  -  \  <  <p  <  \  of  the  axis  0  =  0 .  The  problem  of 


FIG.  1.  The  bubble,  some  stream  lines  of  the  flow,  and  the 
x  and  y  axes  are  sketched.  The  x '  and  /  axes,  used  in  Sec.  V, 
are  also  Indicated. 
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t>,-P<f/2=Pi-ok,  on -i<<p<{,  0=0.  (3) 

Here,  />,,  p,  and  q  are,  respectively,  the  stagnation 
pressure ,  density,  and  the  speed  of  the  fluid  outside 
the  bubble.  In  dimensionless  variables  (3)  becomes 

(?  =  k-y,  on  -  5  <  <?  <  i ,  ,  (4) 

where  the  dimensionless  parameter  y  is  defined  by 

y  »2^-^)/(2<ra)»/V/*  .  (5) 

The  problem  can  be  further  simplified  by  requiring 
the  bubble  to  be  symmetric  about  the  line  y  =x.  This 
implies  that 

y,(-<p,o)=-x,(<p,o) ,  .  (6) 

By  using  Eq.  (6)  we  can  restrict  our  analysis  to  the  in¬ 
terval  0  <<?<£. 


III.  NUMERICAL  PROCEDURE 

Before  solving  the  problem  we  note  that  y¥(<p,  0)  must 
be  singular  at  <p  =  £  like  (i  -  (p)~l/ 2.  Therefore,  we  eli¬ 
minate  this  singularity  by  replacing  <p  with  the  new 
variable  0  defined  by 

<p=32-0*  .  (12) 

Then,  we  introduce  the  AT  mesh  points  0/  given  by 
3,  =  (J-l)/2l/l(N-l),  /=  1,  ...  ,N  ,  (13) 

and  the  N  -  1  corresponding  unknowns 

w . "->•  <»> 

If  we  define  y*v  by  (14)  with  J=N,  it  then  follows  from 
(11)  and  (12)  that 


II.  REFORMULATION 

It  is  convenient  to  reformulate  the  boundary  value 
problem  as  an  integro-differential  equation  by  consid¬ 
ering  the  function 

(<p+  iip)1/2(x  ,  +  iy  ,)  -  (b  /2)1/2 , 

which  is  analytic  in  the  half-plane  $>0  and  vanishes  at 
infinity  as  a  consequence  of  Eq.  (1).  Therefore,  on 
=  0,  its  real  part  is  the  Hilbert  transform  of  its  imag¬ 
inary  part.  By  symmetry  its  imaginary  part  vanishes  on 
^  =  0,  \<p  |>  £  and  therefore  the  Hilbert  transform  yields 


o)-(|Y'M  f1/a 

*  \2/  w  J0  <p’-<p 


i  r 0  (~_r )1/2x^',o) 

"•'-1/2  <P'~<P 

(7) 


We  now  use  the  symmetry  condition  (6)  to  rewrite  (7) 
in  the  form 


/l>\  1/2  ,/>-l/2 

^,(^,0)  =  (-)  (pml/2  +  — — 


:  f  (<p*)l/2yA<p*  f0)  / — - — 4 — - — )  dip* . 
J0  *  \<Pf-<P  <P'  +  <P  J 

(8) 


(15) 


Now,  we  define  0/+i/2:=i(0/  +  0/.1)/2,  and  then  consider 
x6(fiul/z,0).  We  compute  it  from  (8)  rewritten  in  terms 
of  0.  However,  the  integrand  contains  the  factor  (<p'Yf2 
=  (2  -  0*)1/2 ,  which  has  a  singular  derivative  at  0  =  2'1/a. 
In  order  to  integrate  this  factor  accurately,  we  replace 
y9  by  Cy*  -3^  )  +  y*,  use  the  trapezoidal  rule  with  mesh 
point  0/  on  the  intervals  containing  y0  -yj  ,  and  analy¬ 
tically  evaluate  the  integral  containing  y Then, 
x0(0/+!/2,O)  is  given  in  terms  of  they^. 

Next,  we  compute  yfl(0/#iy!l,O),  *M(0Al/2,O),  and 
y  sfl(0/i/2i°)  in  terms  of  the  y’j  by  using  four- point  dif¬ 
ference  formulae.  From  these  expressions  we  calculate 
y,,*„  and  y90  at  (0/+l/2,O),  by  noting  that 


(*+iy)„# 


_(*  +  *y)aa 

~  40s 


Then,  we  substitute  these  expressions  into  (9)  at  the 
N  -  1  points  (0/*i/j,O),  /-l,  .  .  . ,  N-  1  to  obtain  N-l 
nonlinear  algebraic  equations  in  the  JV  + 1  unknowns  yj, 

1  =  1 . N  and  b.  The  last  two  equations  are  obtained 

from  (10)  and  (15)  by  using  three-point  Lagrange  extra¬ 
polation  formulae  to  evaluate  their  left  sides. 


Next,  we  express  the  boundary  condition  (4)  in  terms  of 
x9  and  y,,  noting  that  q2  =b2(x\  +  ya )“l.  Then,  (4)  be¬ 
comes 


b2 

s.+yl 


Lx,.  -  xmymm 


\<P  |<2,  0  =  0  . 


(9) 


Now,  (8)  and  (9)  together  constitute  a  nonlinear  in¬ 
tegro-differential  equation  for  y,(</>)  in  the  interval 
i/>  =  0.  The  symmetry  conditions 

x.ti90)  =  0  ,  (10) 


*,(0,0)*-y.(0,0)  ,  (ID 

complete  the  formulation  of  the  problem  for  yv(<P,0) 
and  b.  This  formulation  of  the  problem,  and  the  num¬ 
erical  method  to  be  described  follows  closely  the  work 
of  Vanden-Broeck  and  Keller.1,2 


The  N  nonlinear  equations  were  solved  by  Newton's 
iteration  method.  For  some  large  value  of  y,  the  initial 
approximation  for  the  bubble  surface  was  taken  to  be 
a  circular  arc.  Iterations  were  continued  until  the 
solution  converged  within  a  given  tolerance.  This  solu¬ 
tion  was  then  used  as  the  initial  guess  for  a  smaller 
value  of?  and  so  on. 

After  a  solution  for  the  y}  and  b  converged  for  some 
value  of  y,  the  profile  40,0),  y(0,O)  was  obtained  by 
integrating  x$  and  y0  .  Typical  profiles  for  various 
values  of  y  are  shown  in  Fig.  2.  For  y  large  the  bubble 
surface  is  close  to  a  circle  of  dimensionless  radius 
y*1.  For  smaller  values  of  y,  the  bubble  surface  is 
close  to  a  square  with  rounded  corners.  Asy  decreases 
further  the  bubble  surface  elongates  in  the  directions 
midway  between  the  axes,  developing  four  horns  or 
spikes. 


M 
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r  =  0 


FIG,  2.  Computed  bubble  profiles  for  various  values  of  y. 
These  profiles  were  found  by  the  method  of  Sec.  ID. 

The  efficiency  of  the  numerical  scheme  was  found  to 
be  limited  by  the  high  curvature  at  the  ends  of  the  spikes. 
Accurate  solutions  for  y  <  -  1.3  could  not  be  computed 
even  with  N  =  50.  Thus,  in  the  last  section  we  shall  con¬ 
sider  an  analytical  approach  to  determine  the  ultimate 
form  of  the  horns  as  y  is  further  decreased. 

From  the  numerical  solution  we  can  calculate  the 
area  A  of  the  bubble  and  the  potential  energy  V.  The 
potential  energy  V  is  equal  to  the  surface  tension  o 
times  the  length  of  the  surface.  In  dimensional  vari¬ 
ables  we  have 

/0//V/3  r  2*1/2 

V  =  \^T)  8J0  MW.Ol+jflM  )t/2dP,  (16) 

(2/t\2/3  /  r  2-*/2 

7?)  8('J0 

+  .  (17) 

We  also  introduce  the  Weber  number  W  defined  by 

W=A1/a(paa/2o)l/3.  (18) 

The  integrals  in  (16)  and  (17)  were  evaluated  by  the 
trapezoidal  rule.  The  results  are  shown  as  functions  of 
y  in  Fig.  3. 

We  see  that  the  potential  energy  V  is  a  monotone  de¬ 
creasing  function  of  y.  However,  the  dimensionless 
area  A,  which  is  equal  to  the  square  of  the  Weber  num¬ 
ber  W2  f  is  not  monotonic  but  has  a  maximum  at  y  --0.2. 
Thus,  the  figure  shows  that  there  is  a  maximum  or 
critical  value  of  W9  above  which  there  is  no  steady  sol¬ 
ution  of  the  kind  considered  here.  This  maximum  value 
is  W9~2.1.  The  existence  of  a  maximum  Weber  number 
has  previously  been  found  for  a  two-dimensional  bubble 
or  drop  in  a  uniform  potential  flow  by  Vanden-Broeck 
and  Keller,3  and  for  a  two  dimensional  drop  in  a  Stokes' 
flow  by  Buckmaster  and  Flaherty3  and  by  Vanden- 
Broeck.4  In  the  three-dimensional  case  the  correspond¬ 
ing  results  were  obtained  by  Moore,5  El  Sawi,0  and 
Miksis  et  al ,7  for  potential  flow  and  by  Rallison  and 
Acrivos*  for  Stokes'  flow. 

There  are  two  solutions  for  each  value  of  the  Weber 
number  in  an  interval  below  the  maximum  value.  For 
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the  branch  on  the  left  of  the  maximum  in  Fig.  3,  the 
potential  energy  V  decreases  as  W  increases,  so  this 
branch  is  probably  unstable.  For  the  branch  on  the 
right  of  the  maximum,  V  increases  as  W  increases,  so 
it  is  probably  stable. 

IV.  SOLUTION  FOR  7  LARGE 

For  y  large  we  seek  the  bubble  surface  as  a  small 
perturbation  of  a  circle  of  radius  y“l.  Thus,  we  write  it 
in  polar  coordinates  r ,  0  as 

r=y“1+y-4p(0)  +  O(y*5)  .  (19) 

The  curvature  k  of  this  surface  is  given  by 

k=y  +y-2[p"(0)+  p(0)]  +  O(y-3)  .  (20) 

Now,  the  unperturbed  complex  potential  <pQ+iil> 0  is  given 
as  a  function  of  z=x  +  iy  by 

<Po  +  *A  =  *2/2+l/2yV  .  (21) 

From  (21)  we  see  that  on  the  circle  of  radius  y’1,  </>„ 

=y“2  cos  20  and,  therefore,  q0  =  -2y~l  sin20. 

In  the  boundary  condition  (4)  we  use  Eq.(20)for  fc  and  the 
preceding  result  for  q0.  The  terms  in  y-3  lead  to 


p"(0)+  p(0)  =4  sina20  .  (22) 

The  solution  p  of  (22)  satisfying  pr(0)  =  p'(r/4)  =  0,  when 
used  in  (19),  yields  the  result 

r=y~1  -  2(1  +  ^  cos40)y“4  +  0(y’5)  .  (23) 

We  next  substitute  (23)  into  the  definitions  of  V  and  A 
and  get 

V  =  (2cft/pa!)t''3  2ir(y'1  -  2y"4)  +  0(y"5)  ,  (24) 

A  =  (2<r / pa2)2 ' 3  7r(r'2-4)--5)  +  0(r-6)  .  (25) 


The  formulae  (24)  and  (25)  agree  with  the  numerical 
results  of  Sec.  m  within  0.5%  for  y  >  4. 


A,  W2 


FIG.  3.  Computed  value?  of  the  bubble  area  A  [in  units  of 
(2<t/(X*2)2/3]  and  potential  energy  V  [In  units  of  (2o4/paJ),/3[ 
as  functions  of  y,  based  upon  Eqs.  (16)  and  (17).  With  this 
choice  of  units  the  lower  curve  Is  also  a  graph  of  W2t  the 
square  of  the  Weber  number. 
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V.  THE  LIMITING  SHAPE  OF  THE  BUBBLE 


Wheny  is  less  than  about  -  1.3,  the  bubble  has  four 
slender  spikes.  They  are  oriented  along  the  lines  y 
-  ±x,  as  shown  in  Fig.  2.  Because  the  spikes  are  so 
slender,  their  shape  can  be  found  approximately  by 
using  the  slender  body  theory  for  bubbles  presented  by 
Vanden-Broeck  and  Keller.2  In  lowest  order,  the  flow 
about  a  symmetric  slender  bubble  is  approximated  by 
the  flow  about  a  rigid  plate  lying  along  the  center  line 
of  the  bubble.  In  the  present  case  the  center  lines  of 
the  four  spikes  consist  of  two  straight  line  segments, 
each  of  some  length  2a,  lying  along  the  lines  y  =  ±x. 

We  introduce  the  coordinates  xr ,  y '  with  axes  along 
these  lines,  and  find  the  potential  b(p{x'  ,y')  of  the  flow 
about  these  plates,  requiring  that  at  infinity  b<p  ~(x* 
-y2)/2.  Evaluating  this  potential  on  the  plate  y'  =  0, 
x*  >  0  we  obtain 

b(p(xr ,  0)  -  (a4  -  x,4Y  /2/2  .  (26) 

By  differentiating  (26)  we  find  that  the  flow  speed  q  on 
the  plate  is 

<7(x',  0)  =  x/3(a4  -  x/4)~1/2,  x’^Q.  (27) 

Before  using  q  to  get  the  bubble  shape,  we  shall  deter¬ 
mine  the  half-length  a  of  the  center  lines.  We  do  so  by 
requiring  the  suction  force  F,  exerted  by  the  flow  on 
the  end  of  a  spike,  to  balance  the  surface  tension  2a. 

As  we  see  in  Ref.  9  [p.  412,  Eq.  (6.5.4)],  F  =  7rpA2/4. 
Here,  A  is  the  coefficient  in  the  expansion  b(p  -Ar1^2 
x  cos0/2  in  terms  of  polar  coordinates  with  their  origin 
at  the  end  of  the  plate.  Upon  setting  F  =2a  and  intro¬ 
ducing  dimensionless  variables  we  obtain 

^=4A.  (28) 

From  (26)  we  find  that  A2  =  a3,  so  (28)  yields 

a  =  (4A)1/3  .  (29) 

We  next  use  (27)  for  q  in  (4)  and  approximate  the  cur¬ 
vature  k  by  -  T7sV(*').  Here,  the  equation  of  the  bubble 
is  y 9  =??(*')  for  x\yf  in  the  first  quadrant  of  the  x,y 
plane,  i.e.,  for  x’  >  yf  >  0.  Then  Eq.  (4)  becomes 

V,'  — .  (30) 

At  the  end  of  the  spike  we  require 

*?(<*)  =  0.  (31) 

In  addition,  by  symmetry,  the  bubble  must  be  normal  to 
the  y  axis  where  it  crosses  that  axis.  In  the  primed 
coordinates  this  axis  is  the  line  xf  =y',  and  this  sym¬ 
metry  condition  becomes 


Vs-*,  (32) 

where  x *  =  rj(x'). 

Upon  integrating  Eq.  (30)  twice  and  using  Eq.  (31),  we 
obtain 

~T  (^_x')  log(a-x') 

-  j-(a  +  *')log(a  +  *')-~  logfo2**'2) 

.  oV,  mlx9  at  ,  .  a*  a2 

+  T*“  7“#*  -a)-I J  +  yT 

+~  logZa+jlaeZa2  .  (33) 


FIG.  4.  The  ultimate  form  of  the  bubble  at  y=-1.8,  when 
opposite  sides  of  each  spike  just  touch  one  another  and  enclose 
a  small  sub-bubble  near  the  tip. 


The  integration  constant  0  in  (33)  is  to  be  found  in  order 
to  make  ij(x')  satisfy  (32).  For  each  value  of  y,  (32) 
can  be  solved  numerically  for  (3  by  iteration.  Then, 

(33)  gives  the  approximate  shape  of  the  bubble. 

By  solving  for  0  in  this  way  for  various  values  of  y , 
and  examining  (33),  we  find  that  opposite  sides  of  each 
spike  touch  each  other  at  y  =  -  1.8.  Then,  each  spike 
contains  a  small  sub-bubble  near  its  tip,  as  shown  in 
Fig.  4.  This  profile  represents  the  ultimate  form  of 
our  family  of  solutions.  For  smaller  values  of  y  oppo¬ 
site  sides  of  each  spike  cross  one  another,  which  is 
physically  inadmissible.  We  could  obtain  physically 
acceptable  solutions  fory<-  1.8  by  allowing  the  pres¬ 
sure  in  the  sub-bubble  to  be  different  from  that  in  the 
main  bubble.  We  shall  not  find  them  because  they  are 
probably  unstable.  It  is  likely  that  the  sub- bubble 
would  become  detached  from  the  main  bubble  if  y  were 
decreased  below  -  1.8. 

Finally,  we  shall  indicate  how  the  preceding  slender 
body  theory  can  be  improved  upon  when  it  is  used  to 
find  a  bubble  with  a  sub-bubble  at  the  tip  of  each  spike. 
We  first  replace  each  sub-bubble  by  a  straight  segment 
of  length  l.  Then,  we  calculate  the  flow  around  the  main 
bubble  terminated  by  these  segments,  and  the  shape  of 
the  main  bubble,  using  the  numerical  method  of  Sec. 

HI.  This  yields  a  two-parameter  family  of  solutions, 
with  parameters  y  and  l .  Then,  for  eachy  we  find  the 
l(y)  for  which  the  forces  balance  at  the  tip  of  each  seg¬ 
ment.  Next,  we  solve  the  equation  analogous  to  (30) 
for  the  function  q(xf)  which  gives  the  surface  of  the 
sub-bubble,  requiring  the  surface  position  and  slope 
to  be  continuous  where  the  main  bubble  meets  the  sub¬ 
bubble.  At  last,  we  determine  the  pressure  in  the  sub¬ 
bubble  by  requiring  that  rj  =  0  at  the  end  of  the  sub-bub¬ 
ble.  If  this  pressure  is  required  to  be  the  same  as  that 
in  the  main  bubble,  we  find  the  value  of  y  corresponding 
to  the  ultimate  configuration  shown  in  Fig.  4. 
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